global oo_ out setup_

%% EXTRACT MATRICES IN DYNARE FORMAT FROM ANDREASEN UNFOLDING
names_dynplus    = vertcat(out.label_y,out.label_v(1:end-M_.exo_nbr,:));  
ny               = size(out.gv,1);
nexo             = M_.exo_nbr;
ns               = size(out.hv,1) - nexo;
nvar             = ny + ns;
select_state     = M_.nstatic+1 : M_.nstatic + M_.npred + M_.nboth;
ys               = oo_.dr.ys(oo_.dr.order_var,:);
setup_.names     = M_.endo_names(oo_.dr.order_var,:); % variables in DR options_.order (needed in DR options_.order since output from policy functions will be in DR options_.order)


%% CONSTANTS
oo_g0_a     = (1/2)*vertcat(out.gss, out.hss(1 : end-nexo,:))...
              + (1/6)*vertcat(out.gsss, out.hsss(1 : end-nexo,:))...
              + (1/24)*vertcat(out.gssss, out.hssss(1 : end-nexo,:));

oo_g0       = 0*oo_g0_a;

for i_index=1:nvar
     name_temp           = setup_.names(i_index,:);
     oo_g0(i_index,:)    = oo_g0_a(strmatch(name_temp,names_dynplus,'exact'),:); %#ok<*MATCH3>
end

oo_.g0_dp = oo_g0;

%% ORDER 1:
oo_g1_a     = vertcat(out.gv, out.hv(1 : end-nexo,:));
oo_g1ssx_a  = (3/6)*vertcat(out.gssv, out.hssv(1 : end-nexo,:));    % only exists if order >2
oo_g1sssx_a = (4/24)*vertcat(out.gsssv, out.hsssv(1 : end-nexo,:)); % only exists if order >3  
       
%% ORDER 2
oo_g2_a     = (1/2)*vertcat(out.gvv, out.hvv(1 : end-nexo,:,:));
oo_g2ssxx_a = (6/24)*vertcat(out.gssvv, out.hssvv(1 : end-nexo,:,:)); % only exists if order >3
       
%% ORDER 3
oo_g3_a     = (1/6)*vertcat(out.gvvv, out.hvvv(1 : end-nexo,:,:,:));

%% ORDER 4
oo_g4_a     = (1/24)*vertcat(out.gvvvv, out.hvvvv(1 : end-nexo,:,:,:,:));
%%

oo_g1       = 0*oo_g1_a;
oo_g1ssx    = 0*oo_g1ssx_a;
oo_g1sssx   = 0*oo_g1sssx_a;
oo_g2       = 0*oo_g2_a;
oo_g2ssxx   = 0*oo_g2ssxx_a;
oo_g3       = 0*oo_g3_a;
oo_g4       = 0*oo_g4_a;


for i_index=1:nvar
     name_temp              = setup_.names(i_index,:);
     oo_g1(i_index,:)       = oo_g1_a(strmatch(name_temp,names_dynplus,'exact'),:); %#ok<*MATCH3>
     oo_g1ssx(i_index,:)    = oo_g1ssx_a(strmatch(name_temp,names_dynplus,'exact'),:); %#ok<*MATCH3>
     oo_g1sssx(i_index,:)   = oo_g1sssx_a(strmatch(name_temp,names_dynplus,'exact'),:); %#ok<*MATCH3>
     oo_g2(i_index,:,:)     = oo_g2_a(strmatch(name_temp,names_dynplus,'exact'),:,:);
     oo_g2ssxx(i_index,:,:) = oo_g2ssxx_a(strmatch(name_temp,names_dynplus,'exact'),:,:);
     oo_g3(i_index,:,:,:)   = oo_g3_a(strmatch(name_temp,names_dynplus,'exact'),:,:,:);
     oo_g4(i_index,:,:,:,:) = oo_g4_a(strmatch(name_temp,names_dynplus,'exact'),:,:,:,:);     
end

oo_.g1_dp     = oo_g1;
oo_.g1ssx_dp  = oo_g1ssx;
oo_.g1sssx_dp = oo_g1sssx;
oo_.g2_dp     = reshape(oo_g2,nvar,(ns+nexo)^2);
oo_.g2ssxx_dp = reshape(oo_g2ssxx,nvar,(ns+nexo)^2);
oo_.g3_dp     = reshape(oo_g3,nvar,(ns+nexo)^3);
oo_.g4_dp     = reshape(oo_g4,nvar,(ns+nexo)^4);

 
%% INCREMENTAL REPRESENTATION (FIRST EXTRACT MATRICES)

%% gxxxx
ghxxxx_t    = vertcat(....
                      out.gvvvv(:    , 1:ns , 1:ns , 1:ns, 1:ns),...
                      out.hvvvv(1:ns , 1:ns , 1:ns,  1:ns, 1:ns));                     
ghxxxx = 0*ghxxxx_t;                     
for i_index=1:M_.endo_nbr
     name_temp                        = setup_.names(i_index,:);
     ghxxxx(i_index,:,:,:,:,:) = ghxxxx_t(strmatch(name_temp,names_dynplus,'exact'),:,:,:,:,:); 
end          
oo_.dr.ghxxxx  = reshape(ghxxxx,nvar,ns^4);
%%
                         
%% guuu
ghuuuu_t = vertcat(...
                   out.gvvvv(:    , ns+1:end , ns+1:end, ns+1:end, ns+1:end),...
                   out.hvvvv(1:ns , ns+1:end , ns+1:end, ns+1:end, ns+1:end));                    
  
ghuuuu = 0*ghuuuu_t;                     

for i_index=1:M_.endo_nbr
    name_temp                = setup_.names(i_index,:);
    ghuuuu(i_index,:,:,:,:,:) = ghuuuu_t(strmatch(name_temp,names_dynplus,'exact'),:,:,:,:,:); 
end      
     
oo_.dr.ghuuuu     = reshape(ghuuuu,nvar,nexo^4);
%%
     
     
%% gxxxu
merge     = vertcat(out.gvvvv,out.hvvvv(1:end-nexo,:,:,:,:));
ghxxxu_t  = merge(:, ns+1:end,1:ns,1:ns,1:ns);     
ghxxxu    = 0*ghxxxu_t;                     

for i_index=1:M_.endo_nbr
    name_temp                = setup_.names(i_index,:);
    ghxxxu(i_index,:,:,:,:,:) = ghxxxu_t(strmatch(name_temp,names_dynplus,'exact'),:,:,:,:,:); 
end      
     
oo_.dr.ghxxxu     = reshape(ghxxxu,nvar,(ns^3)*nexo);
     
     
     
%% gxxuu
merge     = vertcat(out.gvvvv,out.hvvvv(1:end-nexo,:,:,:,:));
ghxxuu_t  = merge(:, ns+1:end, ns+1:end,1:ns,1:ns);     
ghxxuu    = 0*ghxxuu_t;                     

for i_index=1:M_.endo_nbr
    name_temp                = setup_.names(i_index,:);
    ghxxuu(i_index,:,:,:,:,:) = ghxxuu_t(strmatch(name_temp,names_dynplus,'exact'),:,:,:,:,:); 
end      
     
oo_.dr.ghxxuu     = reshape(ghxxuu,nvar,(ns^2)*(nexo^2));
     
   
%% gxxuu
merge     = vertcat(out.gvvvv,out.hvvvv(1:end-nexo,:,:,:,:));
ghxuuu_t  = merge(:, ns+1:end,ns+1:end,ns+1:end,1:ns);     
ghxuuu    = 0*ghxuuu_t;                     

for i_index=1:M_.endo_nbr
    name_temp                = setup_.names(i_index,:);
    ghxuuu(i_index,:,:,:,:,:) = ghxuuu_t(strmatch(name_temp,names_dynplus,'exact'),:,:,:,:,:); 
end      
     
oo_.dr.ghxuuu     = reshape(ghxuuu,nvar,ns*(nexo^3));
%%    


%% ghs2x and ghs2u
merge = vertcat(out.gssv,out.hssv(1:end-nexo,:));
oo_.dr.ghs2x = merge(:,1:ns);
oo_.dr.ghs2u = merge(:,ns+1:end);


%% dr.ghs3x and ghs3u
merge = vertcat(out.gsssv,out.hsssv(1:end-nexo,:));
oo_.dr.ghs3x = merge(:,1:ns);
oo_.dr.ghs3u = merge(:,ns+1:end);
%%


%% ghs2xx, ghs2xu, and ghs2uu
merge     = vertcat(out.gssvv,out.hssvv(1:end-nexo,:,:));
ghs2xx_t  = merge(:,1:ns,1:ns);
ghs2xu_t  = merge(:,ns+1:end,1:ns);
ghs2uu_t  = merge(:,ns+1:end,ns+1:end);

ghs2xx    = 0*ghs2xx_t;
ghs2xu    = 0*ghs2xu_t;
ghs2uu    = 0*ghs2uu_t;

for i_index=1:M_.endo_nbr
    name_temp                = setup_.names(i_index,:);
    ghs2xx(i_index,:,:) = ghs2xx_t(strmatch(name_temp,names_dynplus,'exact'),:,:); 
    ghs2xu(i_index,:,:) = ghs2xu_t(strmatch(name_temp,names_dynplus,'exact'),:,:); 
    ghs2uu(i_index,:,:) = ghs2uu_t(strmatch(name_temp,names_dynplus,'exact'),:,:); 
end    

oo_.dr.ghs2xx     = reshape(ghs2xx,nvar,ns^2);
oo_.dr.ghs2xu     = reshape(ghs2xu,nvar,ns*nexo);
oo_.dr.ghs2uu     = reshape(ghs2uu,nvar,nexo^2);
%%


%% oo_.dr.ghs3 and oo_.dr.ghs4
ghs3_t = vertcat(out.gsss,out.hsss(1:end-nexo));
ghs4_t = vertcat(out.gssss,out.hssss(1:end-nexo));

ghs3    = 0*ghs3_t;
ghs4    = 0*ghs4_t;


for i_index=1:M_.endo_nbr
    name_temp                = setup_.names(i_index,:);
    ghs3(i_index,:,:) = ghs3_t(strmatch(name_temp,names_dynplus,'exact'),:,:); 
    ghs4(i_index,:,:) = ghs4_t(strmatch(name_temp,names_dynplus,'exact'),:,:); 
end 

oo_.dr.ghs3 = ghs3;
oo_.dr.ghs4 = ghs4;
%%

if out.opt.orderApp==1
   oo_.dr.ghxx = zeros(ns+nexo,ns^2);
   oo_.dr.ghuu = zeros(ns+nexo,nexo^2);
   oo_.dr.ghxu = zeros(ns+nexo,ns*nexo);
   oo_.dr.ghs2 = zeros(ns+nexo,1);
end

if out.opt.orderApp>2
    oo_          = full_block_dr_new(oo_,M_,options_);
else
    oo_.dr.ghxxx     = zeros(ns+nexo,ns^3);
    oo_.dr.ghuuu     = zeros(ns+nexo,nexo^3);
    oo_.dr.ghxuu     = zeros(ns+nexo,ns*(nexo^2));
    oo_.dr.ghxxu     = zeros(ns+nexo,(ns^2)*nexo);
    oo_.dr.ghs2x     = zeros(ns+nexo,ns);
    oo_.dr.ghs2u     = zeros(ns+nexo,nexo);
end
